Numerical study of metastable states in the T = RFIM 
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We study numerically the number of single-spin-flip stable states in the T = Random Field Ising 
Model (RFIM) on random regular graphs of connectivity z — 2 and z — 4 and on the cubic lattice. 
The annealed and quenched complexities (i.e. the entropy densities) of the metastable states with 
given magnetization are calculated as a function of the external magnetic field. The results show 
that the appearance of a (disorder-induced) out-of-equilibrium phase transition in the magnetization 
hysteresis loop at low disorder can be ascribed to a change in the distribution of the metastable 
states in the field-magnetization plane. 
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I. INTRODUCTION 

At low temperature, disordered systems may be 
trapped in long-lived metastable states separated by 
large energy barriers that make relaxation to thermal 
equilibrium unlikely on experimental time scales. As a 
consequence, these systems, when driven by an external 
force, display a jerky and hysteretic behavior which is 
associated to irreversible jumps between the metastable 
states, events that are called avalanches. The out-of- 
cquilibrium behavior is thus dictated by the properties of 
the local minima in the energy (or free-energy) landscape 
and it would be useful to know a priori their number and 
energy, or, in the case of magnetic systems, their magneti- 
zation. However, such a computation is a nontrivial task, 
even when the metastable states can be clearly identi- 
fied, as it is the case at zero temperature or in mean-field 
models. Consider for instance the ferromagnetic Random 
Field Ising Model (RFIM) at T = which is a simple pro- 
totype of a spin system with avalanche behavior and has 
been extensively studied in recent years [l[ . In this case, 
one would like to relate the hysteretic response of the sys- 
tem to an externally varying magnetic field to the distri- 
bution of the metastable states in the field-magnetization 
plane. In general, one expects the typical number of 
metastable states to scale exponentially with the system 
size inside the saturation hysteresis loop and to vanish 
outsideQ. At low disorder, however, the existence of 
an out-of-equilibrium phase transition, characterized by 
a macroscopic jump in the hysteresis loop strongly 
suggests that the typical number of states also vanishes 
in some region inside the loop, as depicted schematically 
in Fig. QJb). This feature may also explain the behavior 
observed when the system is driven by the magnetization 
instead of the magnetic field 0, Q. 
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FIG. 1: Putative distribution of the metastable states in the 
field-magnetization plane for the T — RFIM on the cubic 
lattice (or the Bethe lattice with z = 4): (a) strong disorder 
regime (b) low disorder regime. The number of single-spin- 
flip stable states scales exponentially with the system size in 
the shaded region and is zero outside. The solid curve is the 
associated magnetization hysteresis loop. (Color on line.) 



In a recent paperQ, this issue was investigated analyt- 
ically in the case of the RFIM on random regular graphs 
of fixed connectivity z (akin to the Bethe lattice [□]) 
for which the out-of-equilibrium phase transition occurs 
when z > 4 [?[. The annealed complexity T,A{m;H) and 
the quenched complexity T,Q(m;H) of the single-spin- 
flip stable states (i.e., the entropy densities associated to 
the average and typical number of states, Af(m; H) and 
exp[lnA/"(m; H)], respectively) were studied func- 
tion of the magnetization (per spin) m and the external 
field H. The quenched complexity T,Q(m; H), however, 
could only be obtained analytically for z — 2 (i.e., in 
one dimension) and therefore the relationship between 
the phase transition in the hysteresis loop and the dis- 
tribution of the typical states could not be investigated. 
On the other hand, £,4(771; H) was computed for several 
values of z but it was shown that its behavior does not 
reflect the typical behavior of the system; the difference 
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may even be qualitative as for the change in the concavity 
of S J 4(m; if) which already occurs for z = 3. 

In the present work, we revisit the problem by perform- 
ing numerical computations of E ,4(771, if) and Sg(m, if) 
on both random graphs and the cubic lattice. Since ex- 
act enumeration is only possible for systems of very small 
size, we also use more sophisticated statistical methods 
that allow us to consider larger systems and to make sen- 
sible extrapolations of the data to the thermodynamic 
limit despite the fact that the number of mctastable 
states grows exponentially. The case z = 2 for which 
complete analytical results are available is used as a 
benchmark. 

The rest of the paper is organized as follows. In section 
im we introduce the model, define the quantities that we 
want to compute, and briefly recall the main results of 
Ref. 0. In section IIII1 we explain how to properly define 
the complexities in finite systems, and, in section IIV1 
we describe the numerical methods. The results for the 
z = 2 random graphs are presented in section [V] and 
those for z — 4 in section I VII The case of the cubic 
lattice, for which no analytical treatment is possible, is 
studied in section IVlIl Finally, in section fVIIH we give 
some general remarks and conclude. 



II. MODEL 

The RFIM is defined as a set of N Ising spins {s; = 
±1, i — 1, . . . , N} placed on the vertices of a graph, with 
Hamiltonian 

j N N 

U= -- n-tjSiSj - H^Si-^hiSi. (1) 

i 7 j= 1 2—1 i—1 

J is a ferromagnetic coupling constant (we take J = 1 
in the simulations) and n.y is a connectivity matrix that 
defines the topology of the graph (riy = 1 if the vertices 
i and j are connected and otherwise). In the follow- 
ing we shall consider either random regular graphs with 
fixed connectivity z = 2 and z — 4 or the 3d cubic lat- 
tice, if is a uniform external field and the quenched local 
fields hi are independent random variables drawn from 
a Gaussian probability distribution with zero mean and 
standard deviation A which parametrizes the amount of 
disorder in the system. 

Hysteresis and avalanches in the RFIM have been ex- 
tensively studied in the past years using the standard 
zero-temperature Glauber dynamics pi. Ho| . A state 
is then (meta) stable when its energy cannot be decreased 
by flipping a single spin. Although more general dynam- 
ics may be also considered, in which a configuration is 
stable if its energy cannot be decreased by flipping any 
subset of l,2,...,fc spins we shall only focus on the 
1-spin-flip stable states that we simply call metastable 
(all other states are thus unstable with respect to the 
dynamics) 



From Eq. {T]), the change in the energy due to the flip 
of a single spin Sj is 

AH(s t -> - Si ) = 2 Si fi, (2) 

where /j = J^j^i n ij s j + hi + H is the effective local 
field acting on Si. The sum extends over the z spins 
connected with Si (because the random field distribution 
is continuous, the probability of having fi = is zero). 
According to Eq. all the spins Si in a mctastable 

configuration are aligned with their local field, i.e., 

Sifi > 0. (3) 

As will be discussed below, it is also useful to introduce 
a "misalignment" parameter \x (0 < /i < I) defined as 

1 N 

i=l 

where 0(x) is the Heaviside step function (6(cc) = 1 if 
x > and if x < 0). [i is the fraction of spins in a given 
configuration that are not aligned with their local field 
(metastable states thus correspond to fx = 0). 

In the following, the two quantities that describe 
macroscopically a disorder realization are N(m, fi; if), 
the number of states at field if with magnetization 
m = (l/N)J2i s i an d misalignment fi, and J\f(m;H) = 
Jx(m,0;H), the number of metastable states with mag- 
netization m. These quantities are averaged over a 
large number of samples characterized by a different 
set of random fields and random graphs. Since both 
Af(m;H), the average number of metastable states, and 
cxp[ln7V(m; if)], the typical number, are expected to 
scale exponentially with N in some region of the if — m 
plane, the interesting quantities are the corresponding 
annealed and quenched complexities 



Sf(m;ff)= lim -lnM(m;H), (5) 



£g*(m;ff)= lim -In Af(m;H), (6) 

^ N—>oa 1\ 

where the superscript 'oo' indicates that these quantities 
refer to the thermodynamic limit. As shown in Ref. 0, 
the probability of finding metastable states outside the 
average hysteresis loop vanishes when N — > oo. Accord- 
ingly, Eg (m; if) = — oo (or is not defined) in this region 
whereas ET(m;ff) may be positive because of the ex- 
istence of nontypical metastable states. In other words, 
the contour T,^(m; if) = overestimates the size of the 
actual hysteresis loop. The situation inside the loop is 
more complicated, depending on the connectivity z of the 
random graphs (or the spatial dimension) and the disor- 
der strength. For z = 2 (or in one dimension), T,^ and 
Eg are concave functions of if and m for all values of A. 
On the other hand, for z > 3 and A small enough, E^ 3 
becomes a nonmonotonic function of m in some range of 
the field (Fig. 6 in Ref. 0) and m^(H), the (annealed) 



3 



average magnetization of the states (corresponding to the 
maxima of E^(to; H)), may display a discontinuity (Fig. 
9 in Ref . 0) . It is likely that the same is true on euclidean 
lattices in 3 and higher dimensions. The corresponding 
behavior of Eg(m; H) is still unknown but, as discussed 
in Ref. 0, the existence of a jump in the hysteresis loop at 
the coercive field (for z > 4 or d > 3 at low disorder) sug- 
gests that the curve Eg (m; H ) = has a reentrant part, 
as depicted in Fig. [T](b). This implies that Eg(TO;_ff) 
has at least two maxima as a function of to in a certain 
range of the field. It is also possible that the typical mag- 
netization of the metastable states (corresponding to the 
maximum of Eg (m; H)) has a discontinuity at a certain 
value of H . This is the general scenario that we try to 
confirm numerically in the present work. 

III. DEFINITION OF THE COMPLEXITIES IN 
FINITE-SIZE SYSTEMS 

Since numerical simulations arc performed in finite sys- 
tems, one must define the corresponding A^-depcndent 
complexities in such a way that one can properly extrap- 
olate the results to the thermodynamic limit. 

It is quite natural to define the annealed complexity in 
finite systems as 

E A (m; H,N) = ± In Af(m; H, N), (7) 

so that the average number of metastable states is just 
N(m;H 7 N) = exp[iVEU(m; H, N)] . A similar definition 
for the quenched complexity would be Eg (to; 77, N) = 
(l/N)\nAf(m; H, N). This is not convenient, however, 
because the average of the logarithm diverges as soon 
as Af(m; H, N) is zero in some sample, a situation that 
always occurs in a system of small size. One could 
imagine to circumvent the problem by arbitrarily setting 
ln(A/(m; H, N)) = in this case (see, e.g., Ref.[lJ. This 
should be unimportant in the thermodynamic limit since 
it concerns a number of states that is nonexponential in 
system size. However, we recall that we are especially 
interested in the region where the quenched complexity 
is very small and vanishes when N — > oo. This proce- 
dure would thus introduce an unacceptable bias. The 
alternative solution that we propose is inspired by the 
analytical calculation of Ref. |2|, introducing the function 
A a (g; H, N) defined by the following Lcgcndre-Fenchel 
transform 

A a (g;H,N) = max f — In Af a (m;H, N) + gm\ , (8) 
m y N J 

where g € IR (we provisionally add the subscript a to 
recall that Af(m;H,N) is sample-dependent). Note that 
we take the maximum of the expression inside the paren- 
thesis with respect to to, which cures the problem of 
N a {m; H, N) being zero for some values of to. 

Let us first consider the simplest situation where 
Af a (m; H, N) has a single maximum (i.e., is concave). 



E Q (m;H) 
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FIG. 2: Schematic behavior of the quenched complexity in 
the presence of a disorder-induced phase transition (e.g., on 
random graphs with 2 = 4). (a) £<2(m; H) has a single max- 
imum and is positive in the range m_ 00 {H) <m< m +00 (H) 
(b) ~Eq(iti;H) has two maxima and is positive in the range 
mZ^H) < m < ml^H) and m+^H) <m< m+^H). 
m and Sq (m; H) are parametrized by the Legendre variable g 
which increases from — oo to +oo as indicated by the arrows. 

Then, for a given g, the function (l/N) lnAf a {m; H, N) + 
gm (or, equivalently, e N9m J\f a (m; H, N)) has a single 
maximum at some magnetization m a (g; H, N) and 

A a (g;H,N) = — \nJ\f a {m a ,H, N) + gm a . (9) 
The average over disorder yields 



K{g;H,N) = -\nM a {m a ;H,N)+gm (10) 



where m(g;H,N) = m a (g; H, N). This leads to define 
the size-dependent quenched complexity as 

X Q (m;H,N)=A(g;H,N)-gm, (11) 

neglecting terms of order l/N (as one expects that m a 
deviates from it average value by terms of order 1/y/N) 
and considering both to and Eg(m; H, N) as functions of 
the independent parameter g through the relation to = 
m(g; H, N). Clearly, one has Eg (to; H, N) -> Eg°(?n; H) 
when N — > oo, and A and Eg are related by the stan- 
dard Legendre relations dA/dg = to and 3Eg/c?TO = —g 
(neglecting again terms of order l/N to derive the sec- 
ond relation). Varying g from — oo to +oo then traces 
out the curve Eg (to; H, N) vs. m, as depicted schemat- 
ically in Fig. [U^a). Moreover, one readily sees from Eq. 
([HI) that the limit g — ► — oo (rcsp. g — > +oo) selects 
in each sample a the smallest (rcsp. largest) magne- 
tization for which Af a {m; H, N) is strictly positive. As 
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shown in Ref. 0, these two values of to correspond to the 
lower and upper branch of the hysteresis loop, respec- 
tively. Therefore, the two curves obtained by plotting 
m± oc (H, N) = \im g ^± oc m(g; H, N) as a function of H 
define the two branches of the average hysteresis loop. 
On the other hand, the typical (i.e. the most likely) 
value of the magnetization of the metastables states is de- 
scribed by the curve m{g = 0; H, N) which corresponds 
to the locus of the maximum of Eq . 

Now, consider the situation where Af a (m; H, N) has 
two maxima associated to two separated lobes (see e.g. 
Fig. Q2] in section I VI Bp . As discussed previously, such 
a situation is likely to occur in a certain range of H 
when there exists an out-of-equilibrium phase transition, 
for instance on the z = 4 Bethe lattice at low disorder. 
One expects the quenched complexity to display the same 
behavior, as depicted in Fig. 2(b), but this information 
cannot be extracted from Eqs. (fSl fTTj) . In this case, it is 
convenient to treat the two maxima separately by intro- 
ducing two Legendre-Fenchel transforms: 

A~(g;H,N) = max \ —]nJ\f a (m;H,N) + gm 

m<m th y N 

A+{g;H,N) = max ( — lnA/^m; H, N) + gm\ , 

m>m th J 

(12) 

where m t h is a certain threshold in magnetization which 
does not depend on the disorder realization and which is 
chosen so to lie between the two maxima for as many real- 
izations as possible. In consequence, the size-dependent 
quenched complexity also splits into two parts, 



Eq (m;H,N))=A-(g;H,N)-gm- 
E+ (to; H, N)) = A+(g: H, N) - gm A 



(13) 



where to (g;H,N) = m a (g; H, N) and m + (g; H, N) = 

mt(g, H 7 N). It is again useful to consider the mag- 
netization and the complexity Eq or Eq as functions 
of the independent parameter g through the relations 
to = mT (g; H, N) or to = m + (g; H, N). Of course, one 
recovers the preceding formulas with only one function A 
or Eq if | m th |> 1. 

It is also interesting to consider the limits g — > 



±oo. Now m_ oa (H,N) 



lim„ 



to (g; H, N) and 



'+0O 



(H, N) = lim g ^ +00 m + (g; H, N) correspond to the 
lower and upper branches of the hysteresis loop, respec- 
tively. On the other hand, in Eqs. (fl2|) . the limit g — ► +oo 
(resp. g — ► — oo) selects in each sample the largest 
(resp. smallest) magnetization in the interval [— I, mth] 
(resp. [mth,l]) for which Af a (to; H, N) > 0. Therefore, 
if m~(g — > +oo;H,N) < m t h < Wl+(<? — > ~oo; H,N), 
there is a gap in magnetization for which Af a (m; H, N) 
I). On the contrary, if m~(g — * +co; H, N) = m+(g — > 
—oo;H,N), the gap does not exist. In finite systems, 
the width of the gap fluctuates from one sample to an- 
other and the appropriate quantity to be analyzed is the 
average gap. 



IV. NUMERICAL METHODS 

We now detail the numerical methods that we have 
used to compute the number of metastable states as a 
function of the magnetization and the external field. 

A. Exact enumeration 

The simplest method for counting the number of 
metastable states in a given disorder realization con- 
sists in generating all the spin configurations and se- 
lecting those which are metastable in a certain range 
of the field. To do this, we determine the interval of 
stability [H min , H max ] of each configuration s = {si}. 
Eq. (J3j) implies that a spin Si = +1 (resp. s,; = — 1) 
is stable as long as H > —J^2j^TiijSj — hi (resp. 
H < -J^j^riijSj - hi). The fields H min and H max 
correspond to the two situations of marginal stability for 
s, i.e., 



{i|s«=+l} 



{-J^naSj-hi} (14) 



Hmax = min {-J Y^TijjSj - hi} . (15) 
{i|s«=-i} rr • 

3^1- 

Accordingly, configurations with H m i n < H max are sta- 
ble when H is in the range [H m i n , H max ]. On the other 
hand, configurations with H m i n > H max arc always un- 
stable since, whatever the value of H, at least the spin up 
corresponding to H m i n (or the spin down corresponding 
to H max ) is unstable. Since the total number of states 
grows like 2 N , one of course cannot consider large sys- 
tems. In the present work, the largest size investigated 
by this method is N = 26 (with 5000 disorder realiza- 
tions). However, as will be seen below, some of the data 
can be reasonably extrapolated to the thermodynamic 
limit. Moreover, even in a small system, one can observe 
a clear connection between the shape of the hysteresis 
loop and the intervals of stability [H m i n , H max ] of the 
metastable states. 

This is illustrated in Fig. [3] where we plot all the H m in 
and Umax obtained in a single sample (namely, a real- 
ization of the random fields on a random graph of size 
N = 26) for different connectivities z and/or values of A. 
For a given overall magnetization M = Nm, there are in 
general several metastable states, each one characterized 
by a stability range [H m i n , H max ] (for clarity, no link be- 
tween H m in and H max has been drawn in the figure.) The 
comparison of the distribution of these marginal fields in 
the H-m plane with the corresponding saturation hys- 
teresis loop shows that: 

(i) There are no metastable states outside the loop, as 
expected theoretically @- 

(ii) The hysteresis loop is a path in the H-m plane 
that connects the extremal states, i.e. those which have 
the largest H max at a given to (on the ascending branch) 
or the smallest H m in (on the descending branch). An 
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FIG. 3: Limits of stability H min (circles) and H max (crosses) 
of the metastable states in a single disorder realization on a 
random graph of size TV = 26: (a) z = 2, A = 1, (b) z = 2, 
A = 2.5 (c) 2 = 4, A = 1 (d) z = 4, A = 2. 5. (Color on line.) 



avalanche (for instance on the ascending branch) occurs 
between the extremal state s at to and the extremal state 
s' at to' > to when all the stable states with intermedi- 
ate magnetization to < to" < ml have a H max which is 
smaller than H max (s). (In fact, this property is true in 
general for any pair of ordered states s and s' > s 13]). 
In other words, avalanches connect metastable states 
through a sequence of unstable (or, at most, marginally 
stable) intermediate configurations. For z = 2 (and both 
A = 1 and A = 2.5), one can see in Figs. |3[a) and (b) 
that there arc indeed some intermediate marginal states 
along the magnetization discontinuities of the hysteresis 
loop. In contrast, for 2 = 4, marginal intermediate states 
are only observed when A is large enough (Fig.^d)). At 
low disorder (Fig. 02c)), there exists a significant region 
inside the loop where there are no metastable states at 
all. This suggests that the discontinuity in the hysteresis 
loop in the thermodynamic limit is related to the ex- 
istence of a region inside the loop without any typical 
metastable states. 



B. Entropic sampling 

The second method that we have used to compute 
the number of metastable states is the entropic sampling 
Monte Carlo (MC) algorith m [Til , [l5T ] . which is a variant of 
the multicanonical approach [16j. This algorithm explores 
the configurational space with a biased sampling proba- 
bility that favors the weakly degenerate macrostates and 
disfavor the highly degenerate ones. The final result is 
M(m, n; H), with JV(m, (i = 0; H) as a by-product. 

The practical implementation of the method is as fol- 
lows. Starting from an initial guess for Jx(m, fi; H) (for 
instance J\o(m, /i; H ) = 1 for all to and /i), one im- 
proves the estimate iteratively via a sequence of Monte 



Carlo runs. The number of states J\fk(iri, /u; H) ob- 
tained after iteration (or "stage" [l5[) k is used to bias 
the acceptance probability in the next MC run (i.e., 
the transition s — ► s' is accepted with the probability 
min.{l,yvfe(m, n; H) / Afk(m' , //; H)}). One then obtains a 
histogram Hfc + i (to, fi) of the frequency of the macrostates 
and jVfc(m, /i; H) is updated according to the rule 



Nk+i(m, fj,; H) = 



\ Afk(m, n; H) if H fe+ i (m, fi) = 



I A/fc (to, /i] H) ■ Hfc_|_i (to, fx) otherwise. 

(16) 

In principle, this method could yield the actual 
Jx(m, fi; H) after a single very long MC run, even if the 
initial guess is very crude (choosing No{m, fj.; H) = 1 im- 
plies that the states are first sampled at random since 
all spin flips are accepted) . The problem with this brute- 
force procedure, however, is that the number of MC steps 
required to detect the macrostates that are exponentially 
less degenerate than the other ones increases exponen- 
tially with N. In consequence, it is better to update 
N (to, /i] H) iteratively so to disfavor the most degener- 
ate macrostates and facilitate the search of those which 
are less degenerate. We have followed the strategy pro- 
posed in Refs. [IJ, [l5| which consists in starting with 
short MC stages and increasing the length of the subse- 
quent stages if the number of visited states decreases. In 
this way, the number of visited states increases and even- 
tually saturates after a certain number of iterations k* 
that mainly depends on the system size and the length 
of the MC stages. Saturation indicates that all the exist- 
ing macrostates have been already visited (and the corre- 
sponding histogram is flat). In this situation, the number 
A/fe* (to, [i] H) corresponding to non- visited states is set to 
(instead of the initial guess 1) and one finally obtains 
M(m, [i; H) from A4* (m, fi; H) via the rescaling 



A*(m,ti;H), (17) 



N 



which ensures that the sum ^2 m M is equal to 2 

The method is illustrated in Fig. Q] which shows the 
exploration of the states in the m-fj, plane at different 
MC stages. In this case, the number of MC steps during 
the iteration procedure varied from 10 4 initially to typi- 
cally 10 5 during the last ~ 15 stages. By definition, all 
the states are contained in the domain T> = {(to, fj,)\m £ 
[— 1, l],fj, <E [0, 1]}. During the first stages (Fig. 0£a) and 
(b)), only the states close to m = and /i = 0.5 are vis- 
ited. Those close to the boundaries of V are less degen- 
erate and are only visited after a certain number of iter- 
ations (see Fig.[4jc)). Note that certain regions of T> are 
not visited because to and /i are not independent quanti- 
ties and the two coupled equations to = (1/iV) Sj and 
fi = (l/N)J2i®(— s ifi) (with N unknowns {si}) may 
have no solution at all. When there are solutions, which 
correspond to the states visited by the algorithm, their 
number is precisely J\f(rn,/j,;H). Consider for instance 
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the situation at the boundaries of V. For m = ±1, only 
one macrostate is visited (see Fig.[4fc)). Indeed, only the 
value ^ = (1/N) 6(T/i) is compatible with m = ±1. 
Therefore, Af(±l,n,H) = 1 for this value of [i and zero 
otherwise. Imposing fj, = or fj, = 1 is less restrictive 
because fixing fj, only imposes some inequalities on the 
products Sifi (for instance, /i = implies Sj/j > 0, Vi). 
Therefore, in general, several macrostates are visited on 
these two boundaries. 



=5,0.5- 




FIG. 4: Entropic sampling algorithm applied to a single disor- 
der realization on a random graph of size N = 100 with z = 2 
(A = 2 and H = 0). The black dots indicate the states in 
the m-fi space that are visited during the iteration procedure 
after (a) 1, (b) 20, and (c) 60 MC stages. The gray area is 
the domain of existence of all possible states. 



C. Simulated annealing 

As introduced in Sec. IIII1 the magnetization 
m a (g;H,N) which corresponds to the maximum of 
e Ngm J\f a (rn; H, N) for a given disorder realization plays 
a crucial role in our analysis of the distribution of the 
typical mctastable states in the H-m plane. In this sec- 
tion, we describe a Monte Carlo procedure that yields 
m a (g; H, N) without exploring the whole (m, /i) space. 
A similar procedure has been proposed in the context of 
granular materials [ljj and spin systems [IU . 

By sampling the spin configurations with the transi- 
tion probability p(s — > s') = min{l, cxp(Ng(m' — m))}, 
one would obtain a distribution for the magnetizations 
which is proportional to e 9Nm Af(m, //; H, N). According 
to the arguments of the previous section, the maximum 
of this distribution is dominated by the most degenerate 
states which correspond to [i > 0. This simple procedure 
is therefore inappropriate to analyze the statistics of the 
metastable states (/i = 0). In order to sample properly 
these states, one needs to introduce a bias parameter 
such that the transition from s to s' is accepted with 
the probability p(s — ► s') = min{l, exp[N(A(m', fi') — 
A(m 1 /i))]} where A(m,g) = gm — g^y. The interesting 
limit for visiting the metastable states is then g^ — > oo. 
However, since these states are much less degenerate than 
those with /i > 0, taking a large value of g^ does not 
lead to the maximum of e 9Nm Af(m; H, N) right away. 
Instead, one needs to use a simulated annealing algo- 
rithm which consists in performing a sequence of MC 
runs, starting with g^ = and increasing slowly g^ until 
no additional spin flip is accepted (we typically performed 
sequences of 2000 MC runs, increasing g^ by increments 
of 0.05 up to <7 M — 15). 

If J\f a (m; H, N) has a unique maximum, the magneti- 
zation resulting from this procedure is then m a {g; H, N). 
In practice, however, there may be several values of m for 
which e 9Nm Af a (m; H, N) is close to its maximum and the 
final magnetization may slightly differ from the actual 
m a (g; H, N). Our simulations show that such deviations 
do exist but, as will be seen in Sees. IV Bl and IVIBI by 
comparing with the analytical results for z = 2 and the 
results of the other methods, they are not statistically 
relevant. The quenched complexity can then be obtained 
by considering g as a function of m and integrating nu- 
merically the Legcndre equation 9Sq/9?ti = — g, 



The metastable states, which correspond to fx = 0, 
are exponentially less degenerate than the states with 
fi ~ 0.5. Since the algorithm detects the most degener- 
ate states first, the enumeration of the metastable states 
becomes more and more difficult as the system size in- 
creases. However, this method represents a noticeable 
improvement over the exact enumeration procedure and 
we could study sizes up to = 100 (with 300 disorder 
realizations). 



Eq(to) = - / g(m')dm' (18) 

-oo 

using the fact that Eq^-oq) = since there is only one 
state along the hysteresis loop. 

When AT a (m; H, N) has two distinct maxima separated 
by a gap (see e.g. Fig. [12]) so that a threshold m t h has 
to be introduced in the analysis, the method does not 
seem to work. Indeed, we found that convergence to the 
metastable states (expected for g M — > oo) could not be 
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FIG. 5: (a) Annealed complexity T, A (m;H,N) for H = 
and various system sizes (random graphs with z = 2 and 
A = 2). Open symbols: exact enumeration; solid symbols: 
entropic sampling. The dashed line corresponds to the an- 
alytical result of Ref. 2. (b) Dependence of EA.max on the 
system size N. The solid line is a fit to the data. The hor- 
izontal dashed line indicates the exact value EJ maj[ = 0.1 
calculated in Ref. |2j. (Color on line.) 

achieved for all values of g. Nevertheless, it will be shown 
below that some useful informations can be extracted 
from the data even if the two branches m ± (g; H, N) are 
not analyzed separately. With this algorithm, systems 
with sizes up to N = 1000 (with 300 disorder realiza- 
tions) could then be investigated. 

V. RESULTS FOR RANDOM GRAPHS WITH 

2 = 2 

We first present the numerical results for random 
graphs with connectivity z = 2. In this case, the particu- 
lar value of A does not introduce any qualitative change 
and we arbitrarily take A = 2. Our main goal is to check 
the validity of the numerical approaches by comparing to 
the analytical results of Ref. and testing whether the 
different methods are consistent with each other. 



A. Annealed complexity 



for different sizes N. The data for N < 26 have been 
obtained by exact enumeration and those for N > 26 
by entropic sampling. (Recall that the simulated anneal- 
ing method is only appropriate to compute the typical 
quantities.) The good overlap of the two sets of data for 
N = 26 shows that the two methods are consistent. At 
fixed to, the complexity increases with N with a clear 
tendency towards the value S^ 3 (to; H = 0) computed in 
Ref. [H (dashed line in the figure) . Let us analyze in more 
detail the case to = that corresponds to the maximum 
of the complexity at zero field. The size-dependence of 
^U.max plotted in Fig.[5{b) is well described by the ansatz 
£^ max = E^ Ia!iX +aN- 1 lnN+bN- 1 involving the three 
free parameters max , a, and 6[l9|. A least-squares fit 
to the data yields £% max = 0.10, a = -0.46, b = -0.94, 
and the value of E?„„ v is in remarkable agreement with 
the analytical result of Ref. y (actually, a fit of the data 
for N < 26 already provides a good estimate of £j4 max ). 
One can perform a similar analysis of the data for m =/= 
but, because of the discreteness of the overall magneti- 
zation M , an interpolation of the curves £,4(771; H, N) is 
then necessary to study the dependence with N at fixed 

TO. 

We now consider the complexity at H = 1 (see 
Fig. (6]Ja)) to illustrate the performance of our numeri- 
cal procedures in non-zero field. The curves in Fig. [6ja) 
are now asymmetric, with a maximum at a positive mag- 
netization. Note that the entropic sampling method is 
not able to explore the whole range of magnetizations 
in the largest systems. Indeed, when the field H is pos- 
itive, the states with to < are much less degenerate 
than those with to > 0, and are therefore more difficult 
to detect. For instance, using the value of the complex- 
ity for N = 100, the probability of visiting a state with 
to = —0.2 is approximately 5000 times smaller than the 
probability of visiting the states with maximal complex- 
ity (with to ~ 0.42). The detection of rare states is also 
limited by the number of disorder realizations investi- 
gated. This number is only ~ 200 for large systems, 
which is significantly smaller than the 5000 realizations 
generated in the exact enumeration. However, we shall 
show below that the states which are not detected in large 
systems are nontypical and do not affect the quenched 
averages. 

The numerical data for finite systems again underes- 
timate T^{m, H) and the results must be extrapolated 
to N — > 00. The dependence of the maximal complex- 
ity IU,max with N is shown in the inset of Fig. [6] The 
same fit as for H = yields S^ max = 0.079, which is in 
reasonable agreement with the exact value 0.073. 

The curve rriA,max(H) shown in Fig. [5{b) represents 
the locus of the maxima of the annealed complexity in 
the H-m plane. In the thermodynamic limit, the (an- 
nealed) average magnetization of the states defined for 
finite systems as 



We first focus on the annealed complexity as defined by 
Eq. ([7|. Fig.[^a) shows as a function of to at H = 
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FIG. 6: (a) Annealed complexity H, N) for H = 1 and 

various system sizes (random graphs with z = 2 and A = 2). 
The convention for symbols and lines is the same as in Fig. [S] 
The inset shows the size dependence of the complexity at the 
maximum, (b) Average magnetization mA{H) obtained from 
Eq. (|19p (diamonds) and magnetization mA,max(H , N) at the 
maximum of Ea(ti) (stars) for N = 26. The dashed line is 
the analytical result of Ref. |2j and the dotted line represents 
the saturation hysteresis loop calculated in Ref. l7J.(Color on 
line.) 



is exponentially dominated by m^ 3 max (i?) and there- 



fore 



'(H) 



{ (if). In small systems, however, 



rriA(H,N) coincides with inA } max (H, N) at small fields 
only. In particular, m,4^ max (7J, N) reaches saturation be- 
fore rriA(H 7 N). The deviations decrease with system size 
and the two curves are expected to coincide for large N. 
Interestingly, finite-size effects on 171a(H,N) are much 
less important than on TriA,max.(H, N), and rriA{H, N) 
matches the analytical prediction for m°^{H) even for 
small N . 



B. Quenched complexity 

According to the analytical results of Ref. [2, the 
quenched complexity for z = 2 displays a single maxi- 
mum as a function of m. Therefore, the numerical com- 
putations can be based on the definition of Sq given in 
Eq. (jlip . The first test is to check if the different nu- 
merical methods provide a good estimate of m(g; H, N). 
In Fig. we plot m(g; H = 0, N) and m(g; H = 1, N). 
For N = 26, we compare the results obtained by ex- 
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FIG. 7: Magnetization m(g) for H — and H = 1, and var- 
ious system sizes (random graphs with 2 = 2 and A = 2). 
Open, solid, and shaded symbols correspond to the data ob- 
tained from exact enumeration (EE), entropic sampling (ES), 
and simulated annealing (SA), respectively. The dashed lines 
represent the analytical results of Ref. @-( Color on line.) 
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FIG. 8: (a) Quenched complexity Eq(?tj) at H = for various 
system sizes (random graphs with z — 2 and A = 2). The 
large squares correspond to N = 1000 and the other symbols 
and lines refer to the same sizes as in Fig. [5] The dashed line 
is the analytical result for TV — > oo. The inset shows the size 
dependence of EQ imax , the maximal value of the complexity. 
The solid line is a fit to the data. The horizontal dashed line 
indicates the theoretical value Eg max — 0.082 calculated in 
Ref. 0( Color on line.) 



act enumeration and by entropic sampling, and the good 
agreement between the two curves shows that the two 
methods are consistent. The results for = 1000 only 
correspond to the simulated annealing method. The de- 
viations of the data from the analytical result of Ref. [H 
are due to finite-size effects, and the agreement becomes 
very good for N — 1000. Even more important is the fact 
that the curves converge towards the correct limit when 
g — > ±oo (actually, m(g;H, N) is essentially constant 



when 



> 



1) in spite of the lack of data for negative 



magnetizations in large systems (see Fig. [Ha)). Such a 
result confirms that (i) the states with negative m which 
arc not detected by the entropic sampling method in large 
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systems are not typical, and (ii) the range of magnetiza- 
tions for which the density of the typical states is finite 
(i.e., Eq > 0) is well captured by the numerical methods. 

Fig. |5Ja) shows the quenched complexity Eq(to) for 
H = 0, computed from Eqs. (flU)) and (fTTj) . or, for 
N = 1000 (with the simulated annealing method), from 
Eq. (fT8|) . As could be anticipated from the preceding 
figure, the range of magnetizations for which Eg > 
does not depend significantly on N. In contrast, the val- 
ues of Eq(t7i) in this interval display significant finite- 
size effects. As for the annealed case, the actual com- 
plexity in the thermodynamic limit is underestimated 
but the curves extrapolates to the correct limit. For 
instance, the extrapolation of the maximal complexity 
yields Eg max = 0.082 which is in excellent agreement 
with the value computed in Ref. 0. 



VI. RESULTS FOR RANDOM GRAPHS WITH 

z = 4 

In this section devoted to random graphs with z = 4, 
we shall only focus on the features of the metastable 
states in the low disorder regime (A < A c ) for which 
the hysteresis loop exhibits a jump discontinuity, as dis- 
played in Fig. [9] The behavior for A > A c is indeed 
qualitatively similar to the one discussed above for z = 2. 
Specifically, we shall take A = 1 which is smaller than 
the critical value A c = 1.781260- It has been shown 
that the jump in magnetization corresponds mathemati- 
cally to the appearance of a new stable root in the self- 
consistent equation for the fixed-point probability de- 
rived in Ref. 0- This equation has three real solutions 
in the range H\ < H < H2 (hereafter, for simplicity, 
we only consider the lower branch since the upper one 
can be obtained by symmetry) , with two of them merg- 
ing and becoming complex at H\ and i?2- As discussed 
in Ref. 0, it is tempting to associate the reentrant curve 
corresponding to the "unphysical" root of the equation 
in the range H\ < H < H2 to the limit of existence of 
the metastable states, i.e., to the contour Eq(?ti, H) = 0. 

However, before discussing this issue which concerns 
the quenched quantities (the only ones that are actually 
observable), let us first summarize the results of our nu- 
merical analysis for the annealed quantities. This will 
provide an additional comparison with the analytical re- 
sults of Ref. S 



A. Annealed complexity 

The behavior of the annealed complexity in zero field 
is displayed in Fig. ITOT a) for different system sizes. It 
is similar to the one observed for z = 2 and the func- 
tion remains concave in the whole range — 1 < m < 1. 
Moreover, the results nicely extrapolates to the thermo- 
dynamic limit. For instance, as shown in the inset of 
Fig. llOf ah E^max, the complexity at the maximum, ex- 
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FIG. 9: Saturation hysteresis loop on random graphs with 
z = 4 and A — 1. The solid line is the curve obtained in 
the thermodynamic limit from the equations of Ref. 0- The 
intermediate "unstable" parts (in the range Hi < H < H2 
for the ascending branch) are included. The hysteresis loop 
has a jump discontinuity at H = ±H2 (dashed lines). The 
symbols represent the loops obtained numerically in 1000 dis- 
order realizations of size N = 100 (for clarity, the points are 
not connected). (Color on line.) 
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FIG. 10: Annealed complexity Z A (m;H,N) for H = (a) 
and H = 1 (b) and various system sizes (random graphs with 
z — 4 and A = 1). The convention for symbols and lines is the 
same as in Fig. [5] The insets show the size dependence of the 
maximal complexity £yt,max and the corresponding fits. (Color 
on line.) 



trapolates to the value SJ4 max = 0.183 which is in excel- 
lent agreement with the theoretical value of Ref. 0- 

The behavior in non-zero fields is qualitatively differ- 
ent. In this case, the complexity is not concave in the 
whole range of magnetizations and two local maxima may 
be present. This is illustrated in Fig. [TOTb) for H = 1. 
In this case, the global maximum of the complexity is lo- 
cated at negative magnetizations although the external 
field is positive. This is the general behavior observed at 
small values of H. When the field is further increased, 
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the peak at positive magnetizations (which is very close 
to to = 1 for this value of A) becomes the global maxi- 
mum. As will be shown below, the quenched complexity 
displays a similar behavior. 

The numerical results for finite systems again nicely 
extrapolate to the exact thermodynamic limit computed 
in Ref. 0. The extrapolated result for XU. max is 0.139 
(inset in Fig. HOf b)). which is in very good agreement 
with the theoretical value S^ max = 0.138. One may note 
that the corresponding magnetization only deviates from 
its limiting value m^ max = —0.41 by amounts smaller 
than 2/N. 



B. Quenched complexity 

Our study of the quenched average is based on the data 
obtained by entropic sampling and simulated annealing 
on random graphs of size ./V = 100 and N = 1000. The 
hysteresis loops obtained in 1000 realizations of the disor- 
der are shown in Fig.[5]and illustrate the typical behavior 
of the magnetization curve in finite systems (note in par- 
ticular that the average value of the field corresponding 
to the jump in to is H = 1.72 whereas the coercive field 
in the thermodynamic limit is H — 1.45). 

Let us first consider the case H = 0. As shown in 
the inset of Fig. [TlT bh the number of metastable states 
Af(m; H = 0, N) has a single maximum at m w 0, which 
implies that one can use a single Legendre transform to 
compute m(g) and Sq(to) (Eqs.([5])-([TT])). As can be 
seen in Fig. 111( a). m(g) increases monotonically from 
to_oo « —1 to m+oo ~ +1 (note the good agreement 
between the entropic sampling and simulated annealing 
methods). These two values, which are extremely close 
(but not equal) to ±1, represent the lower and upper 
bounds of the magnetization of the typical states and 
correspond to the lower and upper branches of the hys- 
teresis loop, respectively (one can see in Fig. [5] that these 
values arc indeed very close to ±1). The correspond- 
ing behavior of the complexity is shown in Fig. ITlTb). 
The curve is of course symmetric with respect to m = 
(one cannot see on the scale of the figure that the slope 
(9Sq(to)/9to is infinite at to = to±oo). 

Although the method based on two distinct Legendre 
transforms is neither appropriate nor necessary in this 
case, one may check that it yields the same results and 
that there is no gap in magnetization for any choice of 
rath in the interval [to.oo, m +oa \. According to the argu- 
ments in Sec. Mil this shows that the density of the typ- 
ical states is strictly positive when to_oo < to < m+oo, 
that is to say inside the hysteresis loop. 

We next consider the case H ^ and focus on the 
interesting region around the jump in magnetization. 
As can be seen in Fig. [T^] that shows the raw data for 
Af(m;H = 1,7V = 100) obtained in 200 samples, it is 
now crucial to introduce two Legendre functions A~(g) 
and A + (g) in order to take into account the possible exis- 
tence of a gap in the magnetization of the typical states, 
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FIG. 11: (a) m(g) and (b) T, Q {m;H,N) for H = (ran- 
dom graphs with z — 4, N — 100, and A = 1). Circles and 
squares correspond respectively to the entropic sampling and 
the simulated annealing method. In (b), the inset displays 
the raw data for N(m; H = 0, N) (obtained with the entropic 
sampling method) in a semi- log scale. (Color on line.) 




FIG. 12: Number of metastable states N(m;H = 1,N) ob- 
tained by entropic sampling in 200 disorder realizations (ran- 
dom graphs with z — 4, N — 100, and A = 1). Note that the 
right lobe is concentrated very close to m = 1. 



which would imply that the curve Eq(to, H) = has a 
reentrant part, as anticipated in Fig. \V[b). Whether or 
not this part of the curve is described by the "unphysical" 
root of the fixed-point equation of Ref. (corresponding 
to the branch between Hi and Hi drawn in Fig. [9|) is an 
independent issue which is discussed below. 

Fig. IT3l summarizes the results for the magnetizations 
to (g) obtained by entropic sampling. In contrast with 
the case H = 0, we find that toT^ is strictly smaller 
than to^^ for the four values of H considered in the fig- 
ure (note that m + (g) is always very close to 1). This 
points towards the existence of a gap and two separate 
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FIG. 13: Magnetizations m~(g) and m + (g) for different ex- 
ternal fields chosen in the vicinity of the jump discontinuity 
(random graphs with z = 4, N = 100, and A = 1). The 
threshold in magnetization is mth = 0.9 in all cases. 



branches m (g) in the thermodynamic limit. Note that 
these curves have been calculated with the threshold 
value mth = 0.9. The results do not depend on m t h, 
however, as long as the value is chosen in an appropri- 
ate interval whose length increases with H. For instance, 
for H = 0.5 and H = 1.3, the approximate intervals are 
[0.85,0.95] and [0.2,0.95], respectively. Outside these in- 
tervals one gets = m^, which simply indicates 
that there is a finite density of typical states with mag- 
netization mth- 

The corresponding complexities Eg and Eg computed 
from Eqs. (fT3|) are shown in Fig. [14] The two-peak 
structure resembles that of the annealed complexity in 
Fig. HUT b). but whereas E^ has a minimum between 
the two maxima, Eg is not even defined in the inter- 
val < m < rrit^. This gap, whose size increases 
as H approaches the coercive field, indicates that there 
is no typical metastable state having a magnetization in 
this range. (Note that Eg and Eg seem to go to zero at 

m = rnTLoQ and m = ml^, respectively. It is possible, 
however, that the complexity is small but finite at these 
borders, in contrast with what happens for m = mZ^ 
and m = m^^, i.e., along the two branches of the hys- 



teresis loop 

In table |U the results for mI M and mZ^ are com- 
pared to the values of the magnetization obtained from 
the solutions of the fixed-point equation of Ref. (the 
other branch Eg(m) is so close to m = +1 that compar- 
ing the two limits rol^ and mt^ with the theoretical 
values is not meaningful). Whereas ml^ is in reason- 
able agreement with the theoretical value for the lower 
branch of the hysteresis loop, JnT^ significantly deviates 
from the solution that describes the "unstable" branch 
drawn in Fig. [5J This may be due, of course, to finitc- 



FIG. 14: Quenched complexity T.Q(m;H,N) for H = 
0.5, 0.75, 1 and 1.3 (random graphs with z = 4, N = 100, and 
A = 1). The two branches correspond to Eg(m) and Eg(m), 
respectively. For H = 0.5, the inset shows a magnification of 
£ + (m). (The threshold in magnetization is m t h = 0.9 in all 
cases.) 



size effects. In particular, the fact that our data pre- 
dicts a discontinuity in the magnetization for H = 0.5 
whereas the fixed-point equation of Ref. has only one 
root at this field is somewhat surprising. However, this 
may also indicate that the "unstable" branch obtained 
analytically does not represent the limit of existence of 
the metastable states in the low disorder regime. As far 
as we know, there is indeed no proof that the two curves 
should coincide. 
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-0.999 


-0.995 ± 0.005 




0.79 ± 0.02 


0.75 


-0.998 


-0.997 ± 0.002 
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0.61 ±0.04 
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-0.994 


-0.991 ± 0.004 


0.33 


0.17 ±0.05 


1.3 


-0.97 


-0.96 ±0.03 


-0.24 


-0.21 ±0.04 



TABLE I: Comparison of mI 00 (_ff,iV) and m+^H, N) with 
the values of the magnetization obtained from the solution of 
the fixed-point equation of Ref. (random graphs with z = 
4, TV = 100 and A = 1). m^ yst and the theoretical 

values corresponding to the ascending branch of the hysteresis 
loop and the "unstable" intermediate branch, respectively (for 
H — 0.5, the unstable solution does not exist). 



In order to study larger systems and check that the 
above conclusions are correct, one has to resort to the 
simulated annealing procedure (this is mandatory in the 
case of the cubic lattice, as discussed below in section 
IVIip . Unfortunately, as noted in section IIVC1 this 
method does not work properly when one introduces 
a threshold m t h so to take into account the fact that 
A/"(m; H, N) has two maxima. In consequence, we use 



12 




-0.6 -0.4 -0.2 0.2 0.4 0.6 
8 



FIG. 15: Magnetization m(g) for H — 1 obtained from a 
single Legendre transform (random graphs with z = 4 and 
A = 1). Circles correspond to the entropic sampling method 
(N — 100) and squares to the simulated annealing method 
(N — 100: large squares, N = 1000: small squares). The 
insets show the normalized distribution of m a (g) obtained by 
simulated annealing for g = 0, 0.05, 0.15 and N = 100. (Color 
on line.) 



one Legendre transform only, which has the unpleasant 
consequence that the two branches m ± (g) are mixed. 
Nevertheless, as shown in Fig. [15] for H = 1, the exis- 
tence of a gap in the magnetization of the metastable 
states has a clear signature in m(g). In order to compare 
with the entropic sampling method, let us first consider 
the case of random graphs of size N = 100. One can 
see that the two methods are in perfect agreement. The 
comparison with the results obtained with two Legendre 
transforms (see Fig. [T5|) shows that m(g) = m~ (g) for 
g < g\ ~ 0.04 and m(g) = m + (g) for g > g 2 — 0.15. In 
the range g\ < g < g 2 , the magnetization varies continu- 
ously but steeply from the branch m~(g) to the branch 
m + {g). This behavior comes from the fact that the dis- 
tribution p(m a ) of the magnetizations m a {g) has a single 
peak for g < g\ (because A~(g) > A+(g) in all disorder 
realizations) or g > g 2 (because A~(g) < A+ (<?)), and 
two peaks for g\ < g < g 2 (because there are some real- 
izations for which A~(<?) < A+(g)). Therefore, in the lat- 
ter case, the "mixed" quantity m(g) does not represent 
properly the distribution of the magnetizations. How- 
ever, the rapid increase of m(g) signals unambiguously 
the passage from one branch to the other one (compare 
for instance with the case z = 2 in Fig. [7]). One can see 
that the slope becomes steeper as one goes from N = 100 
to N = 1000, and, therefore, one expects a discontinuity 
in the thermodynamic limit. We stress, however, that 
this discontinuity has no physical meaning and that two 
branches m (g) coexist from g = — oo to +oo for this 
value of H . 

Finally, in Fig. [16j we show the curve m(g = 0; H) 
which represents the locus of the maximum of the 
quenched complexity in the H — m plane and corresponds 
to the typical magnetization of the metastable states. 
As could be seen already in Fig. 1141 the maximum for 
if < 1.3 is located at negative magnetizations despite 
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FIG. 16: Typical magnetization m(g — 0;H) (solid circles) 
of the metastable states for H > (random graphs with z = 
4, N = 1000, and A = 1). The data are obtained by the 
simulated annealing method. The small squares represent 
the hysteresis loops obtained in 100 disorder realizations (for 
clarity, the points are not connected). The case H < is 
obtained by symmetry. (Color on line.) 



the fact that the external field is positive. This behav- 
ior is somewhat unexpected (see the discussion in section 
5 of Ref. [3) and similar to that of the (annealed) aver- 
age magnetization tua (H) (Fig. 9 in Ref. 0) ■ This is in 
contrast with the behavior observed for z = 2 and for 
z = 4 in the strong disorder regime (although the change 
of behavior does not exactly occur at A = A c ). When H 
is further increased, the maximum of Sq becomes larger 
than the maximum of Eg and m(g = 0; H) jumps to a 
positive value very close to 1 (the arrow in the figure in- 
dicates the approximate average field at which the two 
maxima are equal in systems of size N = 1000). As we 
shall see in the following section, the qualitative behavior 
on the cubic lattice is quite similar. 



VII. RESULTS FOR THE CUBIC LATTICE 

In this section, we present some results for the 
metastable states of the RFIM on the cubic lattice. We 
focus on the low disorder regime where the hysteresis loop 
is discontinuous. Our only objective is to show that there 
exists a gap in magnetization with no metastable states 
for a certain value of H smaller than the coercive field. 
By continuity, this implies that there is a whole region 
inside the hysteresis loop where the number of states is 
zero. Since the entropic sampling method is limited to 
systems which are too small (say, N < 5 x 5 x 5), we 
have used the simulated annealing algorithm and ana- 
lyzed the results with a single Legendre transform, as 
discussed above. The system size is N = 10 x 10 x 10 
with typically 100 to 300 disorder realizations. We have 
chosen A = 1.8 which is smaller than the critical disorder 
A c ~ 2.2 estimated in Refs.@,0|). 
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FIG. 17: Saturation hysteresis loops (small squares) for 100 
disorder realizations on a cubic lattice with TV = 10 3 and 
A = 1.8 (for clarity, the points are not connected). The solid 
circles represent the typical magnetization m(g = 0, H) for 
H > obtained by simulated annealing. (Color on line.) 




FIG. 18: Magnetization m(g) for H — 1 obtained from a 
single Legendre transform (simulated annealing method on 
a cubic lattice with N = 1000 and A = 1.8). The insets 
display the normalized distribution of m a for g — 0.05, 0.08 
and 0.15. (Color on line.) 



In Fig. [T71 wc show the hysteresis loops obtained in 
100 disorder realizations. All of them display a big jump 
in magnetization at a coercive field whose average value 
is H = 1.95. For H = 1, the behavior of the magnetiza- 
tion m(g) is quite similar to the one observed on random 
graphs for z = 4. As shown in Fig. [TSl there is a steep 
increase in m(g) around g « 0.08 (to be contrasted with 
the Id case shown in Fig. [7] or the situation for A > A c ), 
with a displacement of the peak in p(m a ), the distribu- 
tion of the magnetizations m a (g; H, N), from m ps —0.25 
to m k 1. The interval in which two peaks are present 
simultaneously is very small. This suggests that there 
are in fact two branches m ± (g), each one corresponding 
to a distinct branch E + (m) or E _ (m) of the quenched 
complexity as in Fig. 1141 The similarity with the scenario 
on random graphs is also confirmed by the behavior of 
m(g = 0; H), the typical magnetization of the metastable 



states. In Fig. [T71 like in Fig. [16J the typical magnetiza- 
tion first decreases as the field increases (revealing that 
the maximum of the complexity occurs at negative mag- 
netizations), and finally, at a field that is smaller than 
the average coercive field, jumps to a value very close to 
+1. 

Note that the fact that the domain of existence of the 
metastable states (in the H — m plane) for the cubic 
lattice is nonconvex is also suggested by the local mean- 
field calculations done in Ref. [j in which the system is 
driven by the magnetization at finite temperature (see 
Fig. 5 in that reference). 



VIII. SUMMARY AND CONCLUSIONS 

We have studied numerically the distribution of the 
1-spin-flip stable states in the RFIM at T = on ran- 
dom regular graphs and on the cubic lattice. Three 
complementary numerical methods (exact enumeration, 
entropic sampling, and simulated annealing) have been 
used to calculate the annealed and quenched complex- 
ities as a function of the magnetization and the field. 
The three methods have been shown to yield consistent 
results when applied simultaneously to the same sample. 
The entropic sampling algorithm allows one to study sig- 
nificantly larger systems than exact enumeration while 
obtaining interesting informations on the whole configu- 
rational spacefill . The simulated annealing algorithm, 
however, is the only one which is adapted to large sys- 
tems. 

In the case of random graphs, comparison with the 
analytical results of Ref. [3 (for £U(m, H) and Eq(to, H) 
when z = 2 and for S J 4(m,if) when z = 4) shows that 
the numerical data can be successfully extrapolated to 
the thermodynamic limit. In particular, for z = 4, the 
concave character of £,4(771) for H = and the non- 
concavity for H ^ are well captured by the numerical 
results. 

The present study allows us to confirm that in the low 
disorder regime (i.e., for A < A c on random graphs with 
z = 4 or on the cubic lattice), there is a whole region 
inside the hysteresis loop without any typical metastable 
states, as illustrated in Fig. QJb). In this case, the 
quenched complexity £g(m; H) is positive in two distinct 
regions separated by a finite gap in magnetization (for H 
close but smaller than the coercive field). The existence 
of the two branches is also confirmed by the discontin- 
uous behavior of the typical magnetization of the states 
as a function of the field. In other words, a macroscopic 
discontinuity along the hysteresis loop merely indicates 
that there is a "hole" in the distribution of the metastable 
states (i.e. the domain of existence of the metastable 
states in the H — m plane is nonconvex). It is likely that 
this explanation is valid for all driven systems at T = 
(at least when the hysteresis loop is the outer bound of 
all metastable states, which is the case when the coupling 
interactions are ferromagnetic). 
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FIG. 19: Schematic behavior of the quenched complexity 
£Q(m) in the A-H plane for different values of A and H 
in the case where a disorder-induced phase transition exists. 
The solid line represents the coercive field H coer (A) at which 
the macroscopic discontinuity in the magnetization occurs. 
This line ends at the critical point (CP). Points A to D cor- 
respond to the low disorder regime (A < A c ) and points E 
and F to the strong-disorder regime (A > A c ). The diagram 
for H < is obtained by symmetry. 

The relation between the phase diagram of the out-of- 
cquilibrium transition in the A — H plane and the behav- 
ior of the quenched complexity Sg(m;iY) as a function 
of m is summarized schematically in Fig. [19] for different 
values of H (for simplicity, we only consider H > 0) . The 
solid line represents the coercive field H coer (A) at which 
the macroscopic discontinuity in the magnetization oc- 
curs. This line ends at the critical point (CP). In the 
strong-disorder regime (points E and F), the complexity 
has a single maximum that moves towards positive mag- 
netizations as H increases. This behavior is qualitatively 
similar to the one observed in 1 dimension (e.g., on ran- 
dom graphs with z = 2). In the low disorder regime, the 
behavior as H increases is more complicated. Whereas 
Sq(to) has a single maximum at m = when H = 
(point A), there are two maxima when H approaches the 
coercive field (point B). Moreover, the global maximum 



is the one with the smallest magnetization, which is neg- 
ative (although the field is positive). At a certain field, 
close to but smaller than H coer , the maximum on the pos- 
itive side of the magnetizations becomes the global one 
(which thus corresponds to a discontinuity in the typical 
magnetization of the states). At H = H coer (point C), 
all the states have a positive magnetization except one. 
When H > H coer (point D), the complexity has a single 
maximum at a large positive magnetization. 

Note that we have not detailed the passage from one 
maximum to two maxima in Eg(m) as H increases. This 
point is still unclear and deserves further investigation. 
For simplicity, we have also assumed that £<g(m) goes 
continuously to zero when m approaches the borders of 
the gap (i.e., when m — > or m — > m^, see Figs. [5] 

and [T9| . However, this is not mandatory. Note also that 
the behavior for larger connectivities may be more com- 
plicated than the one depicted in Fig. [T9l Indeed, when 
z —* 00, one expects to recover the mean- field behav- 
ior with its distinct S shape (which in turn implies that 
Eq(to; H = 0) must have three maxima when z is large). 
This is in fact a very interesting issue which justifies an 
analytical study of the quenched complexity on random 
graphs in the large- z limit. This study should also permit 
to understand the actual status of the "unstable" branch 
that can be computed from the fixed-point equation of 
Ref. 0- The present numerical study seems to indicate 
that this branch does not coincide with the boundary 
Sq(to; H) = but finite-size effects are too important to 
yield a definite conclusion. 

Finally, it would be interesting (but probably highly 
nontrivial) to study the behavior of the quenched com- 
plexity as the critical point CP is approached. It is in- 
deed remarkable that the essential phenomenology of the 
out-of-equilibrium phase transition (including the critical 
behavior) is encoded in £g(m, H). 
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